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Abstract 


One  set  of  approximate  electron  fluid  conservation  equations  is 


A  two-dimensional,  electromagnetic,  electron  fluid  computer  code 
for  microwave-induced  surface  breakdown  in  air  is  developed.  This  code 
is  based  on  finite  difference  approximations  to  Maxwell’s  curl  equations 
for  TMom  modes  in  cylindrical  waveguides  and  electron  conservation 
equations  for  continuity,  momentum,  and  energy.  An  implicit,  block- 
iterative,  optionally  stiff,  variable  step  integrator  solves  the  global  set  of 
equations.  A  2.856  GHz,  1  MV/m  amplitude,  radial  electric  field  pulse  in 
a  highly  dispersive  waveguide  impinges  on  two  different  model  field 
emission  structures  with  the  same  emission  threshold  (0.2  MV/m)  and  rate 
(10"/(m3-s)).  One  structure  is  a  perfectly  conducting  flat  end  plate;  the 
other  is  a  perfectly  conducting  axial  stub  of  0.413  cm  radius  by  1.69  cm 
length,  centered  on  a  flat  end  plate.  Pressures  vary  from  50-800  torr.  In 
both  cases,  peak  electric  fields  occur  where  expected,  initial  electron 
generation  is  from  field  emission,  and  later  generation  is  from  avalanching 
in  air.  The  location  of  peak  generation  is  at  the  field  maxima.  The 
amount  of  late  time  density  varies  qualitatively  with  pressure  as  theoreti¬ 
cally  expected. 

Introduction 


Subjects  of  great  importance  in  the  design  and  operation  of  high- 
voltage  electrical  equipment  are  electron  emission  and  surface  flashover. 
As  a  step  toward  further  understanding  of  these  phenomena  in  gas- 
insulated,  high-power  microwave  structures,  we  have  developed  a  two- 
dimensional,  electromagnetic,  electron  fluid  computer  code  in  an  r-z 
cylindrical  coordinate  system  with  azimuthal  symmetry.  This  code  is 
intended  for  simulation  of  microwave-induced  surface  breakdown  in  air- 
insulated  cylindrical  waveguide  structures.  Other  insulating  gases  can  be 
considered  by  the  specification  of  collision  frequencies  for  ionization  and 
momentum  and  energy  transfer  to  neutral  molecules.  Although  the  initial 
computational  tests  involve  field  emission  from  purely  metallic  structures, 
the  code  can  be  adapted  for  emission  from  dielectrics. 
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where  v.  is  the  ionization  collision  frequency,  g  is  the  volumetric  electron 
emission  rate,  m  is  the  electron  mass,  Vm  is  the  electron-neutral  molecule 
momentum  transfer  collision  frequency,  U  is  the  electron  kinetic  energy, 
vu  is  the  electron-neutral  energy  transfer  collision  frequency,  U0  is  the 
average  neutral  energy,  and  Sj  is  the  average  neutral  ionization  potential. 
For  air,  U0  is  0.025  eV  and  6;  is  1 4  eV.  Among  the  terms  neglected  are 
magnetic  forces,  atomic  and  molecular  electron  losses,  emitted  electron 
momentum  and  energy,  convected  energy  and  momentum,  pressure 
gradients,  and  heat  flows.  A  Maxwellian  velocity  distribution  is  assumed. 
The  electron  emission  rate  is  zero  except  in  computational  cells  next  to 
emitting  surfaces,  where  it  may  be  non-zero  under  prescribed  conditions. 
The  collision  frequencies  used  for  air  are  those  previously  used  for  one¬ 
dimensional  microwave  propagation  in  the  atmosphere  [1,2]  and  two- 
dimensional  propagation  in  waveguides  [3,4], 


Computational  Geometry 


Governing  Equations 

Maxwell’s  curl  equations  for  the  TMom  mode  electromagnetic 
components  in  a  nonmagnetic,  gaseous  medium  in  a  cylindrical  symmetric 
coordinate  system  in  MKS  units  are 
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Figure  1  shows  the  geometry  for  the  computer  code.  The  radial 
coordinate  extends  from  0  to  rm;  the  axial  coordinate  extends  from  0  to  zra. 
Perfectly  conducting  boundary  conditions,  for  which  the  transverse 


where  FI^  is  the  azimuthal  magnetic  field,  Er  and  Ez  are  the  radial  and  axial 
electric  fields,  |l0  and  eo  are  the  free  space  permeability  and  permittivity,  n 
is  the  electron  density  ,  e  is  the  electron  charge  (-1.6  x  10-19  C),  and  vr 
and  v,  are  the  radial  and  axial  electron  velocities. 


X  Ez,  n,  vp  u 

Figure  1.  Computational  geometry. 
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electric  field  vanishes,  are  imposed  at  the  surfaces  rm  and  zm.  The  compu¬ 
tational  space  thus  defined  represents  the  interior  of  a  cylindrical 
waveguide  or  can  with  the  left  end  open.  The  center  line  of  the  can  is  the 
z  axis.  A  cylindrical,  axially  centered,  perfectly  conducting  stub  of  radius 
rs  and  length  ls  may  join  the  flat  end  plate  at  the  right.  Such  a  stub  is 
shown  in  Fig.  1.  Other  concentric  structures  at  the  end  plate  are  also 
possible. 

A  regular  mesh  of  uniform  axial  spacing  Az  and  radial  spacing  Ar 
partitions  the  computational  space.  The  seven  solution  variables — H^,  Er, 
Ez,  n,  vr,  v  ,  and  U — are  located  to  create  a  staggered  grid.  Because  of  its 
convenience  for  central  differencing  of  Maxwell’s  curl  equations,  this  type 
of  grid  appears  extensively  in  explicit,  finite  difference,  electromagnetic 
scattering  codes,  such  that  of  Yee  [5].  As  shown  by  the  symbols  in  Fig.  1, 
occurs  at  grid  points  with  even  axial  indices  and  odd  radial  indices.  Er 
occurs  at  grid  points  with  both  odd  axial  and  radial  indices.  Ez  exists  at 
points  with  both  even  axial  and  radial  indices.  The  four  electron  fluid 
quantities  occur  at  the  same  points  as  Ez. 

Additional  Boundary  Conditions 

For  the  considered  test  cases,  an  input  radial  electric  field  of  1  MV/m 
amplitude  for  the  TM01  mode  is  specified  as  a  time  function  along  the  left 
grid  boundary.  No  other  condition  is  imposed  there.  Since  the  left 
boundary  has  no  radiation  condition,  a  spurious  reflection  occurs  there 
when  a  reflected  wave  from  the  right  side  impinges  on  this  boundary.  The 
calculation  near  the  left  boundary  is  thus  only  good  for  ~2 x,  where  x  is 
the  wave  transit  time  from  the  left  boundary  to  the  right.  This  defect  can 
be  easily  remedied.  The  upper  boundary  is  presently  assumed  to  emit  no 
electrons.  The  right  boundary  is  perfectly  conducting.  When  an  axial  stub 
is  included,  the  condition  that  all  variables,  except  the  electron  density, 
vanish  is  imposed  within  the  stub.  Within  the  stub,  the  electron  density  is 
a  small  number,  typically  10~3m~3,  constant  in  time.  The  face  and  the 
lateral  surface  of  the  stub  are  perfectly  conducting.  The  face  of  the  stub 
and  the  end  plate  emit  electrons  at  the  constant  volumetric  rate  g  when  Ez 
is  positive  and  equal  to  or  greater  than  some  constant  threshold  value  E[h. 
Due  to  azimuthal  symmetry,  is  zero  along  the  bottom  boundary. 

Numerical  Solution  bv  Finite  Differences 

In  the  mesh  interior,  central  finite  differences  replace  the  spatial 
gradients  in  Eqs.  (l)-(4).  At  the  boundaries,  one-sided  differences 
are  used.  The  result  of  this  discretization  is  a  global  set  of  ordinary 
differential  equations  in  time.  The  implicit,  block-iterative,  optionally 
stiff,  variable  step  software  package  GEARBI  then  integrates  the 
system  in  time  [6]. 

GEARBI  has  several  attractive  features  for  electromagnetic  prob¬ 
lems  involving  time  scale  changes  due  to  avalanche  ionization  or  electron 
emission.  The  variable  integration  step  allows  steps  somewhat  greater 
than  the  Courant  stability  limit  when  wave  propagation  is  the  only  time- 
varying  process.  When  ionization  becomes  significant,  the  solver  can  take 
smaller  steps  of  the  order  of  vr1.  The  code  can  thus  use  a  much  smaller 
amount  of  computation  time  than  a  fixed-step  code  set  to  resolve  solution 
variation  for  a  high  ionization  rate.  The  stiff  option  allows  larger  time 
steps  for  stiff  problems,  that  is,  with  disparate  time  scales  of  varying 
amplitude,  than  can  a  code  without  this  capability.  The  block-iterative 
solution  method  reduces  the  memory  requirement  for  a  given  grid  size  as 
compared  to  a  conventional  direct  solver. 

Test  Problems 

Two  test  problems  involving  field  emission  from  perfect  conductors 
are  considered.  The  emission  threshold  is  0.2  MV/m.  The  emission  rate 
is  1011  electrons/(m3~s),  except  as  noted.  These  values  do  not  represent 
any  real  material  surface,  but  serve  as  a  starting  point  for  testing  the  code. 
In  each  problem,  a  1  ns-long  radial  electric  field  pulse  at  2.856  GHz  enters 
the  left,  open  end  of  a  cylindrical  can  of  inner  radius  4.445  cm.  An  end 
plate  terminates  the  can  at  41 .5  cm  to  the  right.  This  geometry  was  chosen 
for  simulation  of  planned  experiments. 

In  one  problem,  the  vertical  end  plate  is  absolutely  flat.  In  the  other, 
a  cylindrical,  perfectly  conducting,  inwardly  protruding,  axially  centered, 
electric  field  enhancing  stub  joins  the  end  plate.  This  stub  is  1.69  cm  long 


and  0.413  cm  in  radius.  The  grid  for  each  problem  has  200  points  in  the 
axial  direction  and  44  in  the  radial  direction. 

Results  for  the  Flat  End  Plate 

As  the  input  pulse  propagates  from  the  open  end  toward  the  end 
plate,  it  experiences  strong  dispersion  as  a  result  of  the  chosen  radius  and 
frequency.  Figure  2a  shows  the  radial  electric  field  at  1 .45  ns,  just  before 
the  leading  edge  of  the  pulse  strikes  the  end  plate.  Figure  2b  shows  the 


(2b) 

Time  =  1.4500  x  10  9  sec 


Figure  2.  Microwave  electric  fields  just  befor  hitting  the  flat  end  plate: 
(a)  radial  field:  ( b )  axial  field. 

axial  electric  field  at  the  same  time.  Instead  of  2.9  cycles  of  a  1  MV/m 
modulated  square  pulse  filling  about  72%  of  the  axial  space,  somewhat 
more  than  3  cycles  of  a  760  kV/m,  dispersed  pulse  fills  the  entire  axial 
space.  The  high-frequency  components  precede  the  lower  frequencies  to 
the  end  plate. 

As  the  wave  initially  impinges  on  the  plate,  an  axial  electrical  field 
positive  half  cycle  builds  up.  When  this  axial  field  reaches  the  emission 
threshold,  emission  begins,  causing  the  electron  density  abruptly  to 
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Figure  3.  Appearance  of  electron  density  due  to  field  emission  at  the  flat 
end  plate:, (a)  axial  field;  (b)  electron  density. 


exceed  the  background  level  of  10“3m“3.  At  an  0.05  ns  output  interval  and 
300  torr,  the  first  electron  density  enhancement  appears  at  1.70  ns  into 
the  calculation.  Figure  3  shows  the  axial  field  and  the  density  at  this  time. 
Since  Ez  first  reaches  the  threshold  at  one  radial  grid  space  above  the  z 
axis,  the  peak  electron  density  occurs  there.  During  the  initial  positive 
Ez  half  cycle,  the  electron  density  at  the  end  plate  expands  radially  out¬ 
ward  from  0.107  cm  to  1.55  cm.  A  maximum  of  4.1948  m~3  occurs  near 
the  axis. 

During  the  next  0.20  ns,  Ez  drops  below  the  threshold,  experiences  a 
negative  half  cycle,  and  then  returns  to  threshold.  During  this  interval 
when  no  emission  occurs,  the  peak  density  drops  to  4. 1 899  m-3.  Since  no 
other  electron  losses  are  allowed,  this  is  probably  evidence  of  convective 
loss.  Within  the  second  positive  half  cycle  of  Ez,  the  peak  density  in¬ 
creases  to  13.453  m-3,  and  the  radial  profile  expansion  reaches  2.79  cm. 
This  expansion,  as  shown  in  Fig.  4,  occurs  primarily  because  the  radial 
span  of  Ez  at  the  end  plate  above  the  threshold  increases  with  time. 

In  the  second  negative  half  cycle  of  Ez,  the  peak  density  decreases  to 


Figure  4.  Electron  density  at  2.35  ns  at  the  flat  end  plate. 


13.434  nr3.  During  the  third  positive  half  cycle,  the  density  profile 
reaches  3.41  cm  radially,  and  the  peak  value  achieves  25.472  nr3.  The 
electron  density  at  the  end  plate  at  2.35  ns,  just  after  the  third  positive  Ez 
half  cycle,  is  shown  in  Fig.  4.  The  three  radial  modulations  in  the  profile 
outward  from  the  axis  are  due  to  additional  emission  during  the  three 
positive  half  cycles  of  Ez.  The  extent  of  the  density  from  each  half  cycle 
is  shown  in  the  figure.  After  the  third  positive  half  cycle,  the  radial 
density  expansion  proceeds  at  a  slower  rate.  The  final  outer  radius  of 
4.03  cm  occurs  at  2.90  ns.  Thereafter  the  peak  density  increases,  and  the 
profile  suffers  slight  changes  in  shape,  but  the  outer  radius  does  not  vary. 

From  displays  like  Fig.  2a,  the  one-way,  initial  transit  time  of  Er 
down  the  axis  is  t  -1 .50  ns.  The  first  reflection  from  the  end  plate  should 
reach  the  left  boundary  at  2x  -3.0  ns.  A  contaminating,  nonphysical 
return  reflection  from  the  unterminated  left  boundary  should  reach  the  end 
plate  at  -4.5  ns.  Thus,  the  calculation  at  the  end  plate  should  be  good 
until  around  4.5  ns. 

After  2.7  ns,  the  peak  density  also  increases  throughout  negative  half 
cycles  of  Ez  when  the  peak  Ez  is  of  the  order  of  1 .5  MV/m  and  the  half 
cycles  are  longer  as  a  result  of  the  dispersive  delay  of  the  lower-frequency 
pulse  components.  This  density  increase  must  be  due  to  avalanche 
ionization  in  the  low-pressure  air  near  the  axis  around  the  negative  peaks 
of  Ez. 

These  calculations  have  been  repeated  from  50-800  torr.  Table  1 
gives  the  peak  electron  density  at  four  times  for  ten  pressures.  Those 
values  denoted  with  an  asterisk  occurred  near'  the  axis  at  input  of  the  can. 
At  100  and  50  torr,  avalanche  ionization  in  the  air  becomes  increasingly 
significant. 

Table  1 .  Peak  electron  densities  forflie  flat  end  plate. 


Pressure 

Peak  Electron  Density  (nr3) 

(torr) 

1.70  ns 

1.75  ns 

2.35  ns 

3.65  ns 

50 

*9.0232  x  103 

*9.4609  x  103 

*5.5600  x  104 

*5.1500  x  105 

100 

*0.59506 

4.2183 

25.981 

1.3467  x  103 

200 

0.36778 

4.1558 

25.543 

169.89 

300 

0.36765 

4.1948 

25.472 

108.05 

400 

0.36676 

4.1796 

25.473 

97.713 

500 

0.36780 

4.1348 

25.348 

90.402 

600 

0.36690 

4.1463 

25.431 

88.876 

700 

0.36739 

4.1945 

25.504 

88.301 

760 

0.36733 

4.1789 

25.552 

88.125 

800 

0.36767 

4.1993 

25.520 

87.649 
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Table  2.  Peak  Electron  Densities  for  the  Axial  Stub 


(5a) 


Time  =  1.7500  x  10  9  sec 


Figure  5.  Appearance  of  electron  density  due  to  field  emission  at  the  axial 
stub:  (a)  Axial  field ;  (b)  electron  density. 


Results  for  the  Axial  Stub 


For  this  problem  at  300  torr  and  an  output  interval  of  0.05  ns,  the 
first  appearance  of  significant  electron  density  is  at  1.75  ns  at  the  stub 
face.  The  axial  electric  field  at  1.75  ns  is  shown  in  Fig.  5a,  and  the  density 
is  shown  in  Fig.  5b.  The  electron  density  is  greatest  near  the  outer  edge  of 
the  stub,  where  the  enhancement  of  Ez  is  greatest.  The  electron  density  at 
the  stub  face  increases  monotonically  in  time.  This  indicates  that  ava¬ 
lanche  ionization  from  field  enhancement  at  the  stub  occurs  during 
negative  half  cycles  of  Ez. 

Electron  density  first  appears  at  the  end  plate  at  1.95  ns,  which  is 
during  the  second  positive  half  cycle  of  Ez.  This  density  extends  from 
1. 14  to  2.58  cm  in  a  radial  profile  similar  to  a  positive  half  cosine. 

Because  of  the  stronger  axial  field  at  the  face  of  the  stub  and  its 
longer  time  above  the  threshold,  the  electron  density  increases  more 
rapidly  at  the  stub  face  than  along  the  end  plate.  At  the  later  times,  the 


Pressure 

(torr) 

Peak  Electron  Density  (m-3) 

1.75  ns 

2.35  ns 

3.65  ns 

50 

*9.3457  x  103 

*5.4777  x  104 

>4.1303  x  1019a 

100 

2.1630 

196.52 

1.3689  x  1017 

200 

2.1631 

54.934 

2.5735  x  108 

300 

2.1633 

40.510 

1.9417  x  105 

400 

2.1631 

36.168 

7.7705  x  103 

500 

2.1631 

34.649 

1.4557  x  103 

600 

2.1633 

34.195 

555.20 

700 

2.1639 

33.652 

309.36 

760 

2.1636 

33.670 

242.85 

800 

2.1634 

33.774 

213.03 

Table  3. 

Peak  Electron  Densities  for  the  Axial  Stub  Versus 

Emission  Rate  at  700  torr 

Emission 

Peak  Electron  Density  (nr3) 

rate  (m3-s) 

-1 

1.75  ns 

2.35  ns 

3.65  ns 

10" 

2.1639 

33.652 

309.36 

1012 

21.630 

339.77 

3.1067  x  103 

1013 

216.28 

3.3737  x  103 

3.0928  x  104 

density  at  the  stub  face  greatly  dwarfs  that  along  the  end  plate. 

Calculations  with  the  stub  have  been  done  from  50-800  torr.  Table  2 
gives  the  peak  electron  density  versus  pressure  for  three  times.  Once 
again  values  marked  with  an  asterisk  indicate  values  at  the  left  side  of  the 
grid  near  the  z  axis.  The  value  in  the  3.65  ns  column  marked  with  the  “a” 
occurred  at  3.27  ns.  Beyond  this  time,  the  code  slowed  down  to  very  tiny 
time  steps  and  could  not  make  any  progress  in  a  practical  interval  of 
computation  time. 

At  1.75  ns,  shortly  after  emission  begins,  there  is  practically  no 
difference  in  the  peak  electron  density  at  the  stub  face  from  100-800  torr. 
At  2.35  ns,  there  is  very  little  difference  in  the  peak  density  from  300- 
800  torr.  However  by  3.65  ns,  the  difference  at  various  pressures  amounts 
to  many  orders  of  magnitude,  which  is  indicative  of  electron  generation  by 
avalanche  ionization.  As  the  pressure  decreases,  the  threshold  electric 
field  for  breakdown  or  strong  avalanching  decreases,  and  the  ionization 
rate  increases.  Decrease  of  the  pressure  corresponds  to  movement  down 
the  right  side  of  the  Paschen  curve  toward  the  minimum  breakdown 
threshold  at  several  ton-.  The  avalanche-generated  density  is  expected  to 
vary  roughly  exponentially  with  the  product  of  the  average  ionization  rate 
and  the  time  interval.  The  ionization  rate  increases  with  decreasing 
pressure.  Therefore  the  density  at  later  times  should  be  much  greater  at 
lower  pressures  than  at  higher  pressures  for  electric  fields  at  and  beyond 
the  breakdown  threshold. 

Table  3  shows  the  peak  electron  density  for  700  torr  at  three  times 
during  calculations  with  an  increase  in  the  constant  emission  rate.  The 
density  at  each  time  scales  almost  exactly  with  the  emission  rate.  This 
behavior  indicates  that  avalanche  ionization  is  not  important  for  the 
calculational  parameters  at  700  torr.  This  is  not  the  case  at  sufficiently 
lower  pressure. 

Comparison  of  Tables  1  and  2  shows  that  at  late  times  avalanche 
ionization  becomes  important  for  the  flat  end  plate  at  pressures  below 
100-200  torr.  Below  400-500  torr,  ionization  becomes  increasingly 
important  for  the  stub.  This  is  because  the  electric  field  enhancement  of 
the  stub  generates  strong  fields,  which  stay  above  the  emission  threshold 
for  longer  time  intervals  than  for  the  flat  end  plate. 


Conclusions 


Two  initial  test  cases  of  a  two-dimensional,  cylindrically  symmetric, 
electron  fluid  computer  code  for  microwave-induced  surface  breakdown 
in  air  with  model  electron  field  emission  show  that  the  code  behaves 
qualitatively  as  expected  over  a  range  of  decreasing  pressure.  Many  other 
aspects  of  the  code  remain  to  be  tested,  including,  most  importantly, 
careful  validation  against  experimental  data  for  a  number  of  structural 
arrangements.  Thus  far,  the  code  offers  good  promise  for  problems  of 
microwave-induced  surface  breakdown  involving  field  emission  from 
materials  with  gas  insulation. 
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